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Abstract: A new algorithm for calculating the stransverse mass, Mt2, in either sym¬ 
metric or asymmetric situations has been developed which exhibits good stability, high 
precision and quadratic convergence for the majority of the Mt2 parameter space, leading 
to up to a factor of ten increase in speed compared to other Mt2 calculators of comparable 
precision. This document describes and validates the methodology used by the algorithm, 
and provides comparisons both in terms of accuracy and speed with other existing imple¬ 
mentations. 
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1 Introduction 

The avid reader of papers on Mt 2 phenomenology might be aware that a paper with 
a similar-sounding title and abstract was published a few months ago by Lester and 
Nachman[l]. Their motivation (to produce a high precision calculator that does not suffer 
from inaccuracies in certain difficult regions of the Mt 2 parameter space) is a similar one 
to this paper, and the underlying methodologies used share a remarkably similar starting 
point (as both rely on the same key principle^). Where the two papers differ, however, is 
that the key motivation behind the work presented here was to produce as fast an algo¬ 
rithm for evaluating M'j '2 as possible (whilst still of course producing results to a similar 
high level of precision); this has led to quite a different methodology when compared with 
[1] that calculates Mt 2 at least as accurately (if not more accurately - see section 4.1), 
and which uses arguably a more straightforward implementation in terms of the coding 
involved^. 

^It transpires that this author and one of the authors of [1] (CGL) had independently come across similar 
work involving the matrix addition of conics, the key idea behind these two latest implementations. 

^However possibly not as robust an implementation as [1] - see section 3.2. 


- 1 - 





2 Mt 2 and its computation 


The stransverse mass, Mt 2 , [2] is a now well-established kinematic technique for hadron 
collider events that, in its simplest guise, uses the missing transverse momenta of what is 
assumed to be two invisible (for example hypothetical R-parity conserving supersymmet¬ 
ric) daughter particles, along with the transverse momenta of the associated visible parti¬ 
cles/jets, to establish a maximal lower-bound on the mass of the assumed pair-produced 
parent particles that gave rise to the visible and invisible particles. In its original incarna¬ 
tion, the masses of the two invisible daughter particles were assumed to be identical (the 
“symmetric” case), however M'j '2 can be trivially extended to (and is being increasingly 
used in) the “asymmetric” case, where the two daughter masses are assumed to be differ¬ 
ent. Properties and generalisations of Mt 2 have been investigated extensively [3-28] and 
the variable has been used in many experimental searches, especially for supersymmetric 
particles (see for example [29-34]). 

In the general case Mt 2 must be solved numerically^, and although various algorithms 
exist for doing so, the important “ellipse bisection” methodology and associated algorithm 
developed by Cheng and Han in 2008 [37] (who, following the notation of [1], will be referred 
to here as “CH”) is regarded as the simplest and most efficient way to calculate Mt 2 (note 
a useful library of implementations, including the CH algorithm, can be found at [38]). In 
the CH implementation, kinematic constraints bound the “trial” mass values for Mt 2 into 
two (one for each “side” of an event) conic regions (elliptical when the visible daughter 
particles are massive, parabolic when they are massless) in the invisible particles’ phase 
space, and computation of Mt 2 is reduced to finding the smallest trial mass that ensures 
the two conic regions intersect i.e. the smallest trial mass that could physically have given 
rise to both sides of the event. The method for finding the intersection involves bisecting an 
interval, a robust computational method that guarantees an answer to a requested level of 
precision with linear convergence. However, as pointed out in [1], the CH implementation 
suffers from accuracy problems when trying to calculate Mt 2 in the areas of its parameter 
space where events have very light, but nonetheless massive, daughter particles as opposed 
to massless ones; that is cusp regions of parameter space where the conic boundaries 
crossover from elliptical to parabolic (a parabola being an ellipse in the limit of infinitely 
large axes). It should be pointed out that this is not an error in the CH algorithm as such; 
any implementation for calculating Mt 2 has to find a way to deal with these problematic 
cusp regions, and CH decided the simplest way was to have an arbitrary minimum mass 
below which the event was regarded as massless, thus leading to reduced accuracy for this 
type of event. A detailed discussion of the various problematic cusp regions can be found in 
Walker’s excellent tour de force on Mt 2 computation [27]; Walker’s paper also provides an 
all-encompassing algorithm for calculating M'j '2 reliably in a variety of “difficult” situations, 
however as noted by [I], due to the ensuing complexity of the algorithm it is considerably 

®In some special cases analytic expressions do exist; for example in the “fully massless case”, where the 
visible and invisible daughter particles are all massless, it has been shown [35] that Mt 2 can be calculated 
from the roots of a quartic equation (see also [36] for further discussions of, and a fast, fully-analytic 
implementation for, Mt 2 in the fully massless case.) 
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slower than the CH implementation. 

The recent release of the algorithm developed by Lester and Nachman detailed in [1] (we 
will refer to this going forward as the “LN” implementation) is an important one. Firstly 
it is designed to be robust in all regions of Mt 2 parameter space; it accomplishes this feat 
using two techniques: the first is that it represents the CH kinematic constraint equation for 
each “side” of an event by a generalised conic equation (technically this is in matrix form); 
this has the advantage of allowing the underlying geometry of the event to move seamlessly 
from elliptical to parabolic (as the masses of the invisible daughter particles approach zero) 
in a smooth manner; the second, and this is the key achievement, is that the algorithm 
does not look for an intersection of the perimeters of the two bound conical regions (as 
per previous implementations), but instead looks for an intersection of the bound regions 
themselves (what LN refer to as the “interiors” of the conics). This ensures that there 
is only ever one solution, that is when the bound regions hrst touch.Secondly, the LN 
implementation can compute both symmetric and asymmetric events whereas, for purely 
historical reasons, the CH implementation can only be used in symmetric situations.^ The 
LN implementation is still a bisection method (and thus still linear in convergence) and 
so it should allow Mj '2 to be calculated for every physical event to a high precision. Even 
better, the LN implementation is about two times faster than the CH method to the typical 
precision of the CH algorithm (~ 0.002 GeV) and only at worst about 50% slower when 
calculating Mt 2 to a precision of r\j 10-12 QgY 

3 Route to a faster algorithm 

3.1 What is meant by a “fast” algorithm? 

As the CH and LN implementations use a bisection method, both algorithms converge 
linearly and thus the computational time is proportional to the precision required (LN 
refer to a proportional relationship between computational time, r, and precision, D, in 
an algorithm i.e. r oc H, as a “fast” method).® What is proposed here is an algorithm 
which exhibits mainly quadratic, but at worst superlinear, convergence in finding Mt 2 to 
a similarly high precision as the LN implementation; this greater convergence rate gives a 
computational speed which is np to ten times faster than any known existing methods for 
evaluating Mt 2 - 

^Note this also allows so-called “unbalanced” situations, where one ellipse is entirely bounded by the other 
(and which in existing implementations requires a separate computational check), to be found automatically 
as in this case the “interiors” of the conics always touch (leading for example to a value of Mt 2 of exactly 
zero in the fully massless case). 

^Indeed some less than successful attempts to modify CH’s algorithm to deal with asymmetric events 
(see Walker’s Table 1 in [27]) was another key motivating factor behind the LN implementation. 

®Of course the actual computational time depends as much on making as accurate a guess as one can 
of upper and lower limits for Mt 2 in a specific event (and thus giving as small a range to be bisected as 
possible) as it does on the linear convergence of the method. For example a clever addition to the LN 
algorithm, the so-called “deci-section optimisation” feature, gives, according to its authors, roughly a factor 
of three increase in speed for Mt 2 values near an event’s kinematic minimum. 
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3.2 The faster methodology 

As per the work of CH, evaluating M'j '2 is reduced to finding the intersection of two ellipses 
(or parabolae, or combination of the two). To do this we will use a method due to [39]; 
although originally developed for ellipsoid intersection', this was then extended to ellipses 
culminating in the lemmas and key theorem (Theorem 6) derived in [40]; although only 
proven for ellipses, there is no reason to believe the theorem should not be applicable to 
parabolae intersection as well, and so going forward we will refer to the theorem as if it 
was derived for general conics (see Appendix A for a detailed analysis of why the theorem 
should be applicable to parabolae as well as to ellipses). 

We begin by defining the equation of a conic region V for side one of an event (note 
in the following we ignore hyperbolic conic regions as these have no relevance for the com¬ 
putation of Mt 2 )'- 


'P ApX^ Bpxy -\- Cpy'^ + DpX + Epy + T < 0 (3-1) 

Where for convenience the transverse momentum coordinates for side one {pix,Piy) 
are written simply as {x,y). In matrix form this can be rewritten as X^PX < 0, where 
(x, y, 1) and P = [T)j] is the following 3x3 real symmetric matrix 


P = 


/ A Pp_ PP 
■^P 2 2 

Bp ^ Bp 

~T ~T 
\ Pp pp jr 
\ 2 2 


(3.2) 


Note P would represent an elliptical region if det{P 2 , 2 ) > 0 (^ 2,2 being the 2x2 sub¬ 
matrix of P with the 3rd row and 3rd column removed), and it would represent a parabolic 
region if det{P 2 ^ 2 ) = 0 (and would be a hyperbolic region if det{P 2 , 2 ) < 0). Similarly the 
conic region for side two of the event in matrix form is defined as Q : X^QX < 0. 

Given these two conic regions, the cubic polynomial /(A) = det{\P — Q) (where A is 
some arbitrary factor) is called the characteristic polynomial and 


/(A) = det{XP - Q) = 0 (3.3) 

the characteristic equation of V and Q. The key result (Theorem 6) of [40] upon which 
the new implementation presented here (and indeed the implementation of LN) is based is 
as follows: given two conic regions V : X^PX < 0 and Q : X^QX < 0 

1. V and Q touch externally if and only if /(A) = 0 has a negative real double root; 

2. V and Q are separate if and only if /(A) = 0 has two distinct negative real roots. 

It is also the case (see Lemma 1 of [40]) that in both these cases (and the case where 
the conic regions overlap and there are no negative real roots - see figure 1) the third root 
of /(A) = 0 must be positive (and real). 

^Interestingly this method was developed in the context of robotics and virtual reality environments 
where the intersection of ellipsoids is a computationally efficient way to detect the collisions of 3D entities 
(represented as a combination of ellipsoids) with each other and with their surroundings. 
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Figure 1. Two elliptical regions and their characteristic polynomial /(A). Overlapping iff /(A) = 0 
has no negative real root, touching externally iff /(A) = 0 has a double negative real root, and 
separate iff /(A) = 0 has two distinct negative real roots. The single positive root is also evident in 
each case (reproduced from [40]). 

These properties of the characteristic equation /(A) = 0 lead to specific relationships 
between the coefficients of /(A) and it is these relationships which are utilised in the 
LN implementation to find the initial intersection point of two conic regions. The new 
implementation presented here however takes a different approach. 

Before outlining the new approach, it is worth briefly reminding ourselves of some key 
features of CH’s kinematically bound regions: 

1. These regions correspond to kinematically allowed transverse momenta for the two 
invisible particles of an event, taking into account: an assumed parent particle mass, 
the energy and momenta of the associated visible particles/jets, the assumed masses 
of the invisible particles, and the event’s overall missing transverse momentum. As 
the assumed parent mass is increased, the kinematically allowed transverse momenta 
regions increase. M'j '2 is the minimum parent particle mass that satisfies the kine¬ 
matic constraints of both sides of the event, and hence will manifest as an initial 
intersection of the two allowed regions in the invisible transverse momentum coordi¬ 
nate system of one side of an event; 

2. In an elliptical scenario, the allowed region begins as a point at a parent mass equal to 
the sum of the visible and assumed invisible mass (the kinematic minimum) and ex¬ 
pands from this starting point as the assumed parent mass is increased; in a parabolic 
scenario, the allowed region begins as a straight ray (a line going from a fixed point 
going out to infinity in some direction) at the same kinematic minimum parent mass, 
and unfolds on either side of the ray into a parabola as the parent mass increases, 
with the vertex of the parabola at the original fixed point of the ray. 


- 5 - 









The above is an embarrassingly brief description of these conic regions, the reader is 
reminded that an excellent and detailed description can be found in Walker’s paper [27]. 
The point of this brief description is to make clear that, as the conic regions are obviously 
functions of the trial parent mass, the resulting characteristic equation is also a function 
of the trial parent mass (as well as A), this dependence being found in the coefficients of 
the characteristic equation. We will thus redefine the characteristic polynomial as f{X;6), 
where 6 is defined as: 


<5 = 






2El 


(3.4) 


Here is the “trial” parent mass, /tat the invisible mass, and mi and Ei the mass 
and energy of the most massive or, if massless, the most energetic side of the event.® As 
mentioned above, the key situation is where the two conic regions touch externally and 
this is where /(A; 5) = 0 has a multiple negative root. It is a well-established fact that the 
discriminant. A, of a polynomial with multiple roots equals zero. The discriminant of the 
cubic characteristic polynomial is calculated from its coefficients, and thus the discriminant 
is a function of the trial parent mass. Finding Mt 2 then reduces to finding the lowest 
positive 5 value (and hence from Equation (3.4) the lowest /ry value) that ensures A = 0 
i.e. finding the lowest positive root of A = 0. It is worth noting that the discriminant 
can and does have larger positive roots; in fact in the elliptical scenario the discriminant 
is an eight-order polynomial in 5 (which reduces to a quartic polynomial in the parabolic, 
i.e. massless, scenario®). Any larger positive roots correspond to trial masses where the 
boundary of one conic region touches the other’s boundary beyond the initial intersection 
of the exteriors (see figure 2). 

Reducing the computation of Mt 2 to a simple root-finding procedure (though one 
where it is always clear exactly which is the correct root i.e. it is the delta value which 
ensures the characteristic equation has a double negative root) makes sense because for a 
smooth, differentiable function such as the discriminant, root-finding using numerical tech¬ 
niques is a well-established field with many efficient algorithms available. One of the oldest, 
simplest and fastest is the Newton-Raphson method (which has quadratic convergence, see 
for example [41]), another is the Regula Falsi method (which has guaranteed superlinear 
convergence as long as a modified algorithm is chosen^®, see for example [42]). Both are 
classic methods but both have some well-known limitations. An important feature of ei¬ 
ther is to ensure an iteration begins from as near to a root as possible, and fortunately 
there are some simple techniques for setting bounds on Mt 2 (see for example [27] for a 
detailed discussion of upper and lower bounds on Mt 2 )- By finding sensible bounds on the 


®These quantities are the same as those used in the CH implementation and are also used in the imple¬ 
mentation proposed here. Note also that the quantity 5 is equivalent to Walker’s quantity F in the limit of 
his one-step decay topology. 

®Thus reconfirming the result first derived in [35] that Mr 2 in the fully massless case can be found 
analytically by solving a quartic equation 

number of modified methods are available and have been tested for this implementation; however 
given the straight-forward nature of the discriminant function, one of the simpler modifications, called the 
Pegasus method, was found to be just as efficient as more complicated modifications. 
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Figure 2. Depicting a typical situation where the descriminant function, A{S), has three positive 
roots, i5i, 62 and ^ 3 . These correspond to values where the two ellipses intersect: the two ellipses 
initially intersect at (5i corresponding to the lowest positive root of the discriminant; for larger trial 
masses, ny, the ellipses expand and the two ellipses intersect on their boundaries at two further 
locations, S 2 and ^ 3 . 

trial parent mass prior to starting the iteration, Mt 2 can be computed with around ten 
iterations using the Newton-Raphson method for the vast majority of events to a precision 
of ~ 10“^^ GeV. However in some cases (it is difficult to be precise but perhaps around 
1%) the Newton-Raphson iteration either finds one of the higher roots (identihed due to 
the characteristic equation’s lack of a double negative root) or simply doesn’t converge 
particularly quickly; in this situation the algorithm needs to use a slower (bisection-like) 
method to set revised bounds for the trial parent mass before it can re-compute the correct 
root. As the Newton-Raphson method is not guaranteed to find a root between the two 
new bounds (one of its flaws unfortunately), this time the algorithm uses the Regula Falsi 
method (which is guaranteed to find the correct root albeit with only superlinear conver¬ 
gence). This potential need to switch between two iteration methods evidences the fact 
that this new implementation requires certain checks that the LN implementation does 
not; although it has been tested successfully on millions of different types of events it is 
still possible that there could be some unusual event that might cause the new algorithm 
to fail to find the correct Mt 2 value. These events would hopefully be extremely rare at 
worst, but, for this reason, the LN implementation should still be regarded as the most 
robust method available for calculating Mt 2 - 
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Figure 3. The left hand plot shows the kinematic endpoint, ±0.5 GeV, of a ti Mt 2 distribution 
for one million events. The distribution would be expected to stop at exactly 0 (having subtracted 
the top mass) which it does. The right hand plot shows the difference between these Mt 2 values 
calculated with the proposed implementation and those due to LN, demonstrating how well the two 
algorithms agree - for the vast majority of events (> 97%) the agreement is better than ~ 10“^^ GeV. 


4 Algorithm validation 
4.1 Accuracy 

As discussed previously, the CH implementation, whilst a significant step-forward in Mt 2 
computation, does suffer from accuracy issues as well as being limited to symmetric Mt 2 
events, and so the LN implementation is to be regarded as currently not only the most 
robust but also the most efficient algorithm available for the high precision computation 
of Mt 2 - Thus to validate the accuracy of this new proposed implementation, it was felt 
sufficient to only compare computed Mt 2 values to those obtained using the LN imple¬ 
mentation. 

In the first scenario a toy MC that pair-produces on-shell particles was used to gen¬ 
erate millions of truth-level ti events; Mt 2 was then reconstructed from the (visible) 6-jet 
and lepton masses and momenta and the (invisible) neutrino missing momenta. An accu¬ 
rate implementation for computing Mt 2 would expect to produce values that have a sharp 
cut-off at the top mass. The left-hand plot in figure 3 shows the upper kinematic end¬ 
point of the resulting distribution having subtracted the top mass (assumed to be exactly 
173.0 GeV), demonstrating the expected sharp cut-off at zero. The right-hand plot of fig¬ 
ure 3 compares this implementation’s computed Mt 2 values with that of the LN algorithm, 
demonstrating excellent agreement to a high degree of precision. 

A more challenging class of events for the calculation of Mt 2 would be for those on 
the boundary between elliptical and parabolic behaviour (i.e. massless or near-massless 
events). To produce a suitable dataset, the toy MC was used to produce decays from a 
pair of hypothetical 300 GeV sparticles produced around 200 GeV above threshold. Each of 
the visible and invisible daughter particles from the two decays had a 50% chance of being 
massless, and a 50% chance of having a mass uniformly distributed between 0 and 10 GeV, 
independent of the masses of the other daughter particles in that event. Thus a dataset 
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Figure 4. These plots show the difference between the Mt 2 values calculated with the proposed 
implementation and that due to LN for one million (balanced) events which contain a mixture of 
massless, near-massless and massive but small (< 10 GeV) visible and invisible daughter particles 
(the left hand plot uses the same range as the corresponding plot in figure 3 for ease of comparison, 
the right hand plot shows the full distribution). Although demonstrating again how well the two 
algorithms agree - better than ^ 10“^^ GeV for 93% of the events - there is clearly a broader peak 
and larger tail compared with the tt Mt 2 distribution in figure 3. Note also the handful of events 
with a difference of ~ 10“^ GeV. 


of massless, near-massless and (small but) massive events was generated to investigate 
the proposed algorithm’s accuracy in the crossover region between elliptical and parabolic 
behaviour^^. 

Figure 4 shows the new algorithm’s computed Mt 2 values subtracted from those ob¬ 
tained using the LN algorithm for one million events^^ generated as just described. In 
the left hand plot the range shown has been kept the same as that for the related plot in 
figure 3 to more easily compare them. For the vast majority of events (~ 93%) agreement 
between the two implementations is still better than ~ 10“^^ GeV, however the peak is 
slightly broader and there is clearly more of a tail. The slightly wider peak and more 
extended tail is due to the greater “volatility” in the Mt 2 computation in the boundary 
between elliptical and parabolic behaviour; in these cusp regions a relatively smaller change 
to the trial parent mass, //y, can produce a relatively larger change in the size of the conic 
regions, thus limiting how precisely the value of hy representing the first intersection of the 
conic regions can be determined i.e. even at machine precision, it is difficult to iterate any 
closer, and one has to accept a slightly less accurate determination of Mt 2 - The right hand 
plot of figure 4 shows the full distribution, and a few discrepancies of order rsj 10-2 GeV 
are seen. These handful of events with the larger discrepancies were investigated more 
thoroughly (as were the outliers in the right hand plot of figure 3) and their Mt 2 values 
were determined by other means. These events it seems are generally those on the cusp of 
being unbalanced; it was found that the proposed implementation’s computed values were 
always closer to the true value of Mt 2 than those due to LN’s implementation in these 

^^This set of events would thus be deliberately similar to those analysed in Figure 4 of LN’s paper. 

^^Note only balanced events were used in this comparison, unbalanced events being trivial to calculate. 
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Figure 5. These plots show the distribution in the relative increase in speed for the new algorithm 
compared with the LN implementation. The left hand plot is for the tt decay events; the right 
hand plot for the massless and near-massless events on the cusp between elliptical and parabolic 
behaviour. 


“quasi-unbalanced” situations (see Appendix B for a detailed analysis of two of the larger 
outliers). It is worth remembering however that the events in figure 4 were artificially 
created to be “difficult”, and thus even with the slightly broader distribution of agreement, 
both implementations still represent excellent accuracy in the field of M'j '2 computation. 

4.2 Speed 

The key motivation behind the proposed implementation was to find not only a more accu¬ 
rate way to compute Mt 2 but, as importantly, as fast a method as possible. As previously 
mentioned, given that root-hnding algorithms have benefited from intense scrutiny over a 
very long period of time, and should therefore represent some of the most efficient numer¬ 
ical techniques available, the assumption has been that if the calculation of M'j '2 could be 
reduced to a root-finding problem, then this would be a sensible approach to trying the 
achieve the goal of improved speed. 

It is difficult to assert a definitive speed for an algorithm given the myriad of hardware 
and software architectures upon which it might be run; as the LN implementation and 
this new algorithm represent two of the most accurate methods available for computing 
Mt 2 , the results presented here will instead focus on the relative performance of the two 
implementations when calculating Mt 2 to a precision of order r\j 10-12 GeVi3. 

The authors of [1] suggest the LN algorithm when running on a modest laptop has eval¬ 
uation rates of ~ 200kHz when evaluating Mj '2 to a precision of ~ lO-^^ GeV; this general 
speed, applicable to the vast majority of events evaluated using the LN implementation, 
was conhrmed in the analysis performed hereini'i. The left hand plot of figure 5 shows the 
speed improvement for the new algorithm over the LN implementation for 25,000 tt decay 

^®The programs used for this speed comparison were compiled using the GCC g++ compiler, with -03 
optimisations turned on, and were run within the Cygwin linux emulator on a 2.4Ghz Intel Core i5 laptop 
running Windows 7. 

^^Note the LN algorithm was always run with the deci-section optimisation feature enabled. 
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events as previously described (where Mt 2 is evaluated 100,000 times for each event). The 
results show that the new algorithm can calculate Mt 2 in excess of 10 times faster than 
the LN implementation for some types of event, that it is never less than around twice 
as fast, and that the average speed increase is around 6-7 times faster. There are three 
distinct regions in the distribution: the region around 10 times faster are for events which 
find the correct root first time using the Newton-Raphson method after only a handful of 
iterations; the region around 7 times are the “typical” events which find the root first time 
with around 10 iterations; finally the more complex region below about 6 times is mostly 
where the Regula Falsi method is required, and the speed of the algorithm is determined 
by how quickly the revised upper and lower bounds can be found prior to using the Regula 
Falsi method. The right hand plot of figure 5 shows the speed increase using the more 
“difficult” massless and near-massless events on the boundary of elliptical and parabolic 
behaviour (as previously described); here 25,000 events were again evaluated 100,000 times 
each. Not surprisingly there is a slight reduction in relative speed; there are fewer events 
in the 10 times faster region, and more in the Regula Falsi region below 6 times, however 
the average speed is only marginally worse at around 5-6 times faster. Thus the conclusion 
with respect to computational speed for the proposed algorithm is that on a standard PC 
one would expect maximum evaluation rates of ~ 2MHz, and average evaluation rates in 
excess of IMHz. 

Finally, it is worth noting that as the bisection method of the LN implementation 
has linear convergence, there is a corresponding increase in speed when calculating Mt 2 
to a lower precision (the authors quote a speed increase of approximately three when the 
precision required goes from ~ 10“^^ GeV to ~ 10“^GeV). However the implementation 
proposed herein operates with mostly quadratic (at worst superlinear) convergence (essen¬ 
tially meaning that when the algorithm gets near to a root, it very rapidly homes in on 
its true value); this means that significantly higher precision is achieved in the algorithm’s 
last few iterations, and so in contrast to the LN implementation, there is not a significant 
speed advantage where a lower precision is requested. 

5 Conclusions 

A high-precision algorithm for the calculation of Mt 2 for both symmetric and asymmetric 
events has been developed which, as well as demonstrating very good stability throughout 
the Mt 2 parameter space, has been shown to be at least as accurate as currently available 
algorithms and can compute Mt 2 up to ten times faster than any previous implementation, 
with average evaluation rates in excess of ~ IMHz on a modest PC, to a precision of the 
order of ~ 10“^^ GeV (comparable to the most precise existing implementation due to [1]). 
It is hoped that this new implementation will be useful in situations where rapid and/or 
multiple computations of Mt 2 are required. A single C++ header hie containing the full 
algorithm discussed herein is available as part of the hrst arXiv submission of this paper. 
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A Extending the validity of the method to parabolae 

As discussed in the main body of the paper, the implementation proposed herein relies on 
the work summarised in [40] , however the proofs presented in that paper focus specifically 

on ellipses (for reasons explained previously). It would be useful to investigate if there is 
any reason to believe that the five key lemmas leading to Theorem 6 of [40] cannot be 
extended to cover parabolae^®. 

The five lemmas are as follows (note a distinction is made between the boundary of an 
elliptical region (i.e. X^AX = 0) and the interior of an elliptical region (i.e. X"^AX < 0), 
and only the part of lemma 5 relevant to the computation of Mt 2 is discussed): 

Lemma 1 The characteristic polynomial derived from two elliptical regions, /(A) = 0, 
either has three positive roots, one positive and two negative roots, or one positive 
and a pair of complex conjugate roots; 

Lemma 2 If the interiors of two elliptical regions do not intersect, then /(A) = 0 has a 
negative root; 

Lemma 3 If the interiors of two elliptical regions intersect, then any real root of /(A) = 0 
is positive; 

Lemma 4 If two elliptical regions touch externally, then /(A) = 0 has a negative double 
root; 

Lemma 5 If /(A) = 0 has a negative double root, then the elliptical regions V : X"^AX < 
0 and Q : X^BX < 0 touch each other externally. 

Lemma 1 and the relevant part of lemma 5 (the second paragraph for those wishing 
to review it) do not distinguish between the type of conic and so the proofs are applicable 
to any conic (indeed the relevant part of lemma 5 is actually proven for general conics). 
Lemma 4 also does not distinguish between the type of conic (however it does depend on 
lemma 2). Lemma 3 does not technically distinguish between the type of conic, but it 
does require there to be a point on the plane exterior to both conics, which for two finite 
bounded ellipses is self-evident. This is not so clear-cut for parabolae, as two parabolic 

does LN’s implementation, since both derive ultimately from the idea to use properties of the 
characteristic polynomial first proposed by [39]. 

^®The following discussion is not attempting to be a formal proof that Theorem 6 can be extended to 
parabolae. However, when read in conjunction with a knowledge of the proofs contained in [40], it is hoped 
the logical arguments put forward in this appendix are sufficient to give the reader comfort that it would 
be relatively straightforward to formally extended Theorem 6 to cover the intersection of parabolae. 
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Q(0) = A 


Q(i) = B 



Figure 6. Configuration for Lemma 2 (reproduced from [40]). 

regions with their vertices infinitely displaced from each other (e.g. one vertex at —oo 
the other at +oo) could in principle together cover the plane. However for any parabolic 
region describing the allowed transverse momenta in the computation of Mt 2 , at least one 
of these parabolae will, for physical reasons, have a finite vertex. This therefore would 
ensure that there is always a point on the plane exterior to both and thus allows lemma 3 
to be trivially extended to parabolae. In order to extend Theorem 6 to parabolic regions 
we are therefore left with only needing to extend lemma 2 , which is slightly more involved 
and is specifically proven for elliptical disks only. We shall take a brief look at the key 
points involved in the proof and then see how they might be extended to parabolae. 

The situation is as depicted in figure 6 , where we have two elliptical regions that do 
not intersect (note however that they are allowed to touch on their external boundary and 
this would not affect the proof). The substitution X = {fj. — l)/fi (which maps /i G [0,1] to 
A G [— 00 , 0 ]) is made which transforms the characteristic equation /(A) = det{XA — B) = 0 
to g{^) = det{{l — g)A + gB) = 0. A new conic can be made, Q{fi) = (1 — fj,)A + fiB, 
and it is straightforward to see that (5(0) is equivalent to the elliptical region A and (5(1) 
is equivalent to elliptical region H; the proof then shows that Q{g) must also represent an 
elliptical region as det{Q 2 , 2 ) > 0 (given that A and B were defined to be elliptical regions 
so that det{A 2 ^ 2 ) > 0 and det{B 2 ^ 2 ) > 0). The crux of the proof is then to show that there 
must exist a point Xi that is interior to an elliptical region Q(/Ui) (where 0 < /ii < 1 ) but 
which is exterior to both A and B. 

To extend this to parabolae we can proceed in a similar way. However now instead of 
defining A and B to be elliptical regions, we allow them to be either elliptical or parabolic 
regions i.e. by allowing det{A 2 fl) > 0 and det{B 2 ^ 2 ) > 0. This then leads, by the same 
argument, to Q(/u) representing either an elliptical or parabolic region i.e. det{Q 2 , 2 ) > 0 . 
If the conic regions are elliptical we already have a proof, if they are parabolic then we 
have the situation depicted in figure 7. Here it is arguably easier to show that there must 
exist a point Xi that is interior to the parabolic region Q(/ii) but which is exterior to both 
parabolae A and B. For example we can translate and rotate A and B (note such affine 
transformations will not alter the roots of the characteristic equation) so that they are 
separated by the y-axis (or if the two conics touch on their external boundary, the y-axis 
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Figure 7. Configuration for Lemma 2 extended to parabolae. 


can be made tangent to both conics at the intersection point). We then have A centered 
at (—oo,—oo), B centered at (oo,—oo), and define Q{ni) to be the parabola centered at 
(0, —oo). It is clear that there is a point Xi located along the y-axis somewhere between 
the vertex and center of Qipi) (but excluding the tangent point with both parabolae if they 
are touching externally) which is internal to Q(/Ui) but external to A and B. This simple 
argument is enough to allow the proof of lemma 2 to be concluded as per the elliptical case 
and thus extend both lemma 2 and thus Theorem 6 to parabolae. 

Although the above is not a formal proof that these lemmas can be extended to the 
parabolic case (and thus that Theorem 6 is valid for parabolae), it should give significant 
comfort that there do not seem to be any obvious impediments to doing so and that a 
formal extension would be a reasonably trivial affair. In their paper, LN also recognise 
the lack of a formal proof (though one of the authors (CGL) similarly conjectures that it 
should be extendable to parabolae), and the above should hopefully give users of the LN 
implementation greater comfort that this conjecture is likely to be true. 

B Analysis of discrepancies in the evaluation of Mt 2 between the Lester- 
Nachman and the proposed implementations 

In section 4.1 we saw in figure 4 that the LN implementation and the one proposed here 
did not always agree on the Mt 2 value, and there were a handful of discrepancies as large 
as ~ 10“^ GeV (and some smaller discrepancies of order ~ 10“® GeV in figure 3). A large 
number of the events which give rise to discrepancies between the two algorithms have 
been looked at in detail and, as they fall into two main “quasi-unbalanced” categories, in 
this appendix we will look in more detail at just two events (chosen to represent the two 
categories); these events are two of the larger differences in Mt 2 value visible in figure 4, 
with discrepancies of 0.00954 GeV and —0.00254 GeV respectively. 

The first event has the following parameters (notation is as before with subscripts 1 
and 2 indicating sides of an event): mi = 0.00473856926 GeV, pix = 110.43799970914 GeV, 
Ply = -213.46687262192 GeV, m 2 = 0.0 GeV, p 2 x = -20.28455035002 GeV, P 2 y = 
235.76522546534 GeV, = 111.30684472357 GeV, = 37.47049084405 GeV, pm = 
PN 2 = 0.0 GeV; the correct Mt 2 value for this event is 0.007098115 GeV. The parameters 
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Figure 8. Discriminant functions for the two events, both showing a zoomed-in section near the 
origin where the lowest root can be found. For the first event the lowest root (5 = 2.4174985 x 
10“^°GeV) equates to a “trial” parent mass, fiy, value of 0.007098115 GeV; for the second event 
the lowest root {6 = 1.4207550 x 10“^ GeV) equates to a /ry value of 0.214304220 GeV. 



Figure 9. Conic regions for the two events. The zoomed-in section shows that the first intersection 
of side I’s ellipse with side 2’s parabola occurs at the trial parent masses obtained from the lowest 
root of the discriminant function, confirming the fj^y values as the correct Mt 2 values for the events. 


for the second event are: mi = 0.21164596081 GeV, pix = 52.78735611764 GeV, piy = 
-34.61581597866 GeV, m 2 = 0.0 GeV, p 2 x = -27.12567045837 GeV, p 2 y = 

- 8.57910177494 GeV, = -245.80036782403 GeV, = -77.06353868081 GeV, pm = 
PN 2 = 0.0 GeV; the correct Mt 2 value for this event is 0.214304220 GeV. 

Figure 8 shows the cubic characteristic polynomial discriminant functions for these 
two events; the first has five real positive roots (only three of which are visible in the plot) 
and the second event has only one real root; the two lowest roots for event 1 (visible in the 
zoomed-in section) are very close together and small in value (i.e. close to zero^^). 

Figure 9 shows the two conical regions for each event, both evaluated at their respective 
trial parent mass valnes as obtained from the lowest root of the discriminant functions 
shown in figure 8. Side 1 of both events is a (very thin) ellipse which first manifests at 

^^Indeed if the small mass of side 1 of this “quasi-unbalanced” event was actually zero, it would be 
unbalanced and the Mt 2 value would be exactly zero. 
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Figure 10. Conic regions for the two events using the Mt 2 values as calculated by the LN 
algorithm. It is clear that event I’s value is too high as side I’s ellipse has grown well beyond the 
point of initial intersection with the parabola of side 2; for event 2 the value is too small as neither 
the ellipse has grown, nor the parabola widened, sufficiently for first intersection. 

the origin, at the respective event’s side 1 mass value. Side 2 for both events is a (very 
thin) parabola whose vertex is the event’s respective missing transverse momenta and 
values^®. It is clear from the zoomed-in sections of the figure that the two conic regions 
do indeed first intersect at the respective trial mass values, thus confirming that these 
values are the correct Mt 2 values for the events. 

For completeness, figure 10 shows the conic regions at the Mt 2 values evaluated by 
the LN algorithm (0.016643073 GeV and 0.211763085 GeV respectively); it is clear that the 
Mt 2 value for event 1 is too large as the ellipse has grown well beyond the point of initial 
intersections with the parabola; for event 2, the value is too low as the parabola needs to 
widen further, and the ellipse expand further, before the conical regions will first intersect. 
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